home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / eulisp / mpfeel.lha / MPFeel / Modules / polly.em < prev    next >
Lisp/Scheme  |  1992-10-06  |  12KB  |  379 lines
  1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  2. ;;                                                                           ;;
  3. ;;  EuLisp Module                     Copyright (C) University of Bath 1991  ;;
  4. ;;                                                                           ;;
  5. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  6.  
  7. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  8. ;                                        ;
  9. ;     C78 Computer Algebra Project I : Polynomial Algebra System        ;
  10. ;                                        ;
  11. ;    Author : Keith Playford ( ma6kjp )                    ;
  12. ;                                        ;
  13. ;    For : Russell Bradford                            ;
  14. ;                                        ;
  15. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  16.  
  17. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  18. ;                                        ;
  19. ;    Polynomial Representation                        ;
  20. ;                                        ;
  21. ;    A sparse recursive form ( analogous to that used within REDUCE )    ;
  22. ;    is used which is canonical with lexicographic ordering or the       ;
  23. ;    variables.                                ;
  24. ;                                        ;
  25. ;    ( See the accompanying notes for a more detailed description )        ;
  26. ;                                        ;
  27. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  28.  
  29. (defmodule polly
  30.  
  31.   (lists
  32.    list-operators
  33.    streams
  34.    extras
  35.    symbols
  36.    strings
  37.    arith
  38.    ccc
  39.    (except (null) class-names)
  40.    generics
  41.    
  42.    threads
  43.    futures) ()
  44.  
  45.   ;; Future arith...
  46.  
  47.   (defmethod binary-plus ((a future-object) b) (binary-plus (futureeval a) b))
  48.   (defmethod binary-plus (b (a future-object)) (binary-plus b (futureeval a)))
  49.  
  50.   (defun mkterm(var exp coef)
  51.     (cond ((null coef) nil)    
  52.       ((zerop exp) coef)    
  53.       (t (cons (cons var exp) coef))))
  54.  
  55.   (defun cons-pol(term f)
  56.     (cond ((null term) f)
  57.       ((numberp term) term)
  58.       (t (cons term f))))
  59.  
  60.   (defun lead-t(f) (car f))
  61.  
  62.   (defun exp-t(term) (cdar term))
  63.  
  64.   (defun exp-l(f) (cdaar f))
  65.  
  66.   (defun var-t(term) (caar term))
  67.  
  68.   (defun var-l (f) (caaar f))
  69.  
  70.   (defun coef-t (term) (cdr term))
  71.  
  72.   (defun coef-l (f) (cdar f))
  73.  
  74.   (defun vex-l(f) (caar f))
  75.  
  76.   (defun remains(f) (cdr f))
  77.  
  78.   ;; Term ordering...
  79.  
  80.   (defun orderp (s1 s2) 
  81.     (or (eq s1 s2) (string-lt (symbol-name s1) (symbol-name s2))))
  82.  
  83.   (defun ordervar(t1 t2)
  84.     (cond ((numberp t1) nil)
  85.       ((numberp t2) t)
  86.       ((not (orderp (var-t t2) (var-t t1))))))
  87.  
  88.   (defun greater-lead-p (vex1 vex2)
  89.     (cond ((eq (car vex1) (car vex2)) (< (cdr vex2) (cdr vex1)))
  90.       (t (not (orderp (car vex2) (car vex1))))))
  91.  
  92. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  93. ;                                        ;
  94. ;    Polynomial Addition                            ;
  95. ;                                        ;
  96. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  97.  
  98.   (defun standadd(a b)
  99.     ((lambda(n) (cond ((zerop n) nil)
  100.               (t n)))
  101.      (+ a b)))
  102.  
  103.   (defun add-pol(f g)
  104.     (cond ((null f) g)
  105.       ((null g) f)
  106.       ((numberp f) (add-con-pol f g))  
  107.       ((numberp g) (add-con-pol g f))  
  108.       ((greater-lead-p (vex-l f) (vex-l g)) 
  109.          (cons (lead-t f) (add-pol (remains f) g)))
  110.       ((greater-lead-p (vex-l g) (vex-l f)) 
  111.          (cons (lead-t g) (add-pol f (remains g))))
  112.       (t (cons-pol (add-term (lead-t f) (lead-t g)) 
  113.                (add-pol (remains f) (remains g))))))
  114.  
  115.   (defun add-con-pol(c pol)
  116.     (cond ((null pol) c)
  117.       ((numberp pol) (standadd c pol))
  118.       (t (cons (lead-t pol) (add-con-pol c (remains pol))))))
  119.  
  120.   (defun add-term(t1 t2)
  121.     (cond ((numberp t1) (standadd t1 t2))
  122.       (t (mkterm (var-t t1) (exp-t t1) 
  123.              (add-pol (coef-t t1) (coef-t t2))))))
  124.  
  125.   (defmethod binary-plus ((a pair) (b pair)) (add-pol a b))
  126.   (defmethod binary-plus ((a pair) (b number)) (add-pol a b))
  127.   (defmethod binary-plus ((a number) (b pair)) (add-pol a b))
  128.  
  129.   (export add-pol)
  130.  
  131. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  132. ;                                        ;
  133. ;    Polynomial Multiplication                        ;
  134. ;                                        ;
  135. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  136.  
  137.   (defun multi-pol(f g)
  138.     (cond ((or (null f) (null g)) nil)
  139.       ((numberp f) (multi-term-pol f g)) 
  140.       ((numberp g) (multi-term-pol g f))
  141.       (t (add-pol (list (multi-term (lead-t f) (lead-t g)))
  142.               (add-pol (multi-term-pol (lead-t f) (remains g))
  143.                    (add-pol (multi-term-pol (lead-t g) (remains f))
  144.                     (multi-pol (remains f) 
  145.                            (remains g))))))))
  146.  
  147.   (defun multi-term-pol(term pol)
  148.     (cond ((null pol) nil)
  149.       ((and (numberp term) (numberp pol)) (* term pol))
  150.       ((numberp pol) (cons-pol (multi-term term pol) nil))
  151.       (t (cons-pol (multi-term term (lead-t pol)) 
  152.                (multi-term-pol term (remains pol))))))
  153.  
  154.   (defun multi-term(t1 t2)
  155.     (cond ((and (numberp t1) (numberp t2)) (* t1 t2))
  156.       ((ordervar t1 t2) (mkterm (var-t t1) (exp-t t1)
  157.                     (multi-term-pol t2 (coef-t t1))))
  158.       ((ordervar t2 t1) (mkterm (var-t t2) (exp-t t2)
  159.                     (multi-term-pol t1 (coef-t t2))))
  160.       (t (mkterm (var-t t1) (+ (exp-t t1) (exp-t t2))
  161.              (multi-pol (coef-t t1) (coef-t t2))))))
  162.  
  163.   (defmethod binary-times ((a pair) (b pair)) (multi-pol a b))
  164.   (defmethod binary-times ((a pair) (b number)) (multi-pol a b))
  165.   (defmethod binary-times ((a number) (b pair)) (multi-pol a b))
  166.  
  167.   (export multi-pol)
  168.  
  169. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  170. ;                                        ;
  171. ;    Polynomial Subtraction                            ;
  172. ;                                        ;
  173. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  174.  
  175.   (defun sub-pol(f g) (add-pol f (multi-pol -1 g)))
  176.  
  177.   (defmethod binary-difference ((a pair) (b pair)) (sub-pol a b))
  178.   (defmethod binary-difference ((a pair) (b number)) (sub-pol a b))
  179.   (defmethod binary-difference ((a number) (b pair)) (sub-pol a b))
  180.  
  181.   (export sub-pol)
  182.  
  183. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  184. ;                                        ;
  185. ;    Polynomial Exponentiation                        ;
  186. ;                                        ;
  187. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  188.  
  189. ; Raises f to the power n ( n a natural number ) using multi-pol by the     ;
  190. ; method of repeated squaring.                            ;
  191.  
  192.   (defun raise-pol(f n) 
  193.     (cond ((< n 0) (error "Attempt to raise poly to a negative power"))
  194.       (t (raise-main f n))))
  195.  
  196.   (defun raise-main(f n)
  197.     (cond ((zerop n) 1)
  198.       ((zerop (remainder n 2)) (raise-main (multi-pol f f) (quotient n 2)))
  199.       (t (multi-pol f (raise-main (multi-pol f f) (quotient n 2))))))
  200.  
  201. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  202. ;                                        ;
  203. ;    Polynomial Pseudo-division                        ;
  204. ;                                        ;
  205. ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
  206.  
  207. ; Tests for unacceptable cases and call the main routine.                   ;
  208.  
  209.   (defun psdiv-pol(u v)
  210.     (cond ((null v) (error "Attempt to pseudo dividefun by zero"))
  211.       ((or (null u) (not (consp v)) (not (consp u))) nil)
  212.       ((neq (var-l u) (var-l v)) nil)
  213.       ((greater-lead-p (vex-l v) (vex-l u)) nil)
  214.       (t (divmain u v (exp-l u)))))
  215.  
  216. ; Recursively constructs the quotient.                        ;
  217. ; Keeps track of the order u should be in case of cancellation.            ;
  218.  
  219.   (defun divmain(u v m)
  220.     (let ((qterm (get-qterm u v m))
  221.       (new-u (get-newu u v m)))
  222.       (cond ((= m (exp-l v)) qterm)
  223.         (t (cons-pol qterm (divmain new-u v (- m 1)))))))
  224.  
  225. ; Returns the pseudo remainder of u v (checking for unacceptable cases)        ;
  226.  
  227.   (defun psrem-pol(u v)
  228.     (cond ((null v) (cerror "Attempt to pseudodivide by zero"))
  229.       ((or (null u) (numberp u) (numberp v)) u)
  230.       ((neq (var-l u) (var-l v)) u)
  231.       ((greater-lead-p (vex-l v) (vex-l u)) u)
  232.       (t (remmain u v (exp-l u)))))
  233.  
  234. ; Works as for pseudo division but does not keep track of the quotient        ;
  235.  
  236.   (defun remmain(u v m)
  237.     (cond ((< m (exp-l v)) u)
  238.       (t (remmain (get-newu u v m) v (- m 1)))))
  239.  
  240. ; Constructs the current quotient term via the given algorithm              ;
  241.  
  242.   (defun get-qterm(u v m)
  243.     (let ((k (- m (exp-l v))))
  244.       (mkterm (var-l v) k (multi-pol (coefn u m (var-l v))
  245.                      (raise-pol (coef-l v) k)))))
  246.  
  247. ; Constructs the new remainder polynomial via the given algorithm           ;
  248.  
  249.   (defun get-newu(u v m)
  250.     (let ((k (- m (exp-l v))))
  251.       (sub-pol (multi-pol (coef-l v) u)
  252.            (multi-pol (multi-term-pol (mkterm (var-l v) k 1) v)
  253.               (coefn u m (var-l v))))))
  254.  
  255. ; Returns the var to the power n th coefficient of f                        ;
  256.  
  257.   (defun coefn(f n var)
  258.     (cond ((null f) nil)
  259.       ((eq var (var-l f)) (cond ((= (exp-l f) n) (coef-l f))
  260.                     (t nil)))
  261.       ((zerop n) f)
  262.       (t nil)))
  263.  
  264. ; The given example polynomials.                        ;
  265.  
  266. (setq r1 '(((x . 5) . 1) ((x . 4) ((y . 1) . 2)) ((x . 2) . 1) ((y . 2) . 1)))
  267. (setq r2 '(((x . 3) ((y . 2) . 1)) ((x . 1) ((y . 1) . 4)) . 1))
  268.  
  269. ; Power product ordering.
  270.  
  271.   (defun power-product-degree (pp)
  272.     (cond ((or (numberp pp) (null pp)) 0)
  273.       (t (+ (exp-l pp) (power-product-degree (coef-l pp))))))
  274.  
  275.   (defun power-product-coef (pp)
  276.     (cond ((null pp) 0)
  277.       ((numberp pp) pp)
  278.       (t (power-product-coef (coef-l pp)))))
  279.  
  280.   (defun power-product-without-coef (pp)
  281.     (cond ((null pp) 0)
  282.       ((numberp pp) 1)
  283.       (t (list (mkterm (var-l pp) (exp-l pp)
  284.                (power-product-without-coef (coef-l pp)))))))
  285.  
  286. ; Need variable dicrimination
  287.  
  288.   (defun power-product-lt (pp1 pp2)
  289.     (< (power-product-degree pp1) (power-product-degree pp2)))
  290.  
  291.   (defun leading-power-product-with-coef (p)
  292.     (cond 
  293.       ((null p) ())
  294.       ((numberp p) 1)
  295.       (t
  296.         (let ((max ()))
  297.       (mapc
  298.        (lambda (term)
  299.          (let ((sub (* (list term)
  300.                    (leading-power-product-with-coef
  301.                      (coef-t term)))))
  302.            (when (power-product-lt max sub) (setq max sub))))
  303.        p)
  304.       max))))
  305.  
  306.   (defun leading-power-product-coef (p)
  307.     (power-product-coef (leading-power-product-with-coef p)))
  308.  
  309.   (defun leading-power-product (p)
  310.     (power-product-without-coef (leading-power-product-with-coef p)))
  311.  
  312.   (defun unity-check (v e c)
  313.     (if (zerop e) c (list (mkterm v e c))))
  314.  
  315.   (defun simple-power-product-quotient (pp1 pp2)
  316.     (cond ((numberp pp1) 1)
  317.       ((numberp pp2) pp1)
  318.       ((eq (var-l pp1) (var-l pp2))
  319.         (unity-check
  320.           (var-l pp1) (- (exp-l pp1) (exp-l pp2))
  321.           (simple-power-product-quotient (coef-l pp1) 
  322.                          (coef-l pp2))))
  323.       (t
  324.         (list (mkterm (var-l pp1) (exp-l pp1)
  325.               (simple-power-product-quotient (coef-l pp1) pp2))))))
  326.  
  327.   (defun power-product-very-simply-divisible-p (pp1 pp2)
  328.     (cond ((and (numberp pp1) (numberp pp2)) (cons pp1 pp2))
  329.       ((or (numberp pp1) (numberp pp2)) ())
  330.       ((and (eq (var-l pp1) (var-l pp2)) (= (exp-l pp1) (exp-l pp2)))
  331.         (power-product-very-simply-divisible-p (coef-l pp1) (coef-l pp2)))
  332.       (t ())))
  333.  
  334.   (defun power-product-lcm (pp1 pp2)
  335.     (cond ((null pp1) pp2)
  336.       ((null pp2) pp1)
  337.       ((and (numberp pp1) (numberp pp2)) (lcm pp1 pp2))
  338.       ((numberp pp2)
  339.         (list (mkterm (var-l pp1) (exp-l pp1)
  340.               (power-product-lcm (coef-l pp1) pp2))))
  341.       ((numberp pp1)
  342.         (list (mkterm (var-l pp2) (exp-l pp2)
  343.               (power-product-lcm (coef-l pp2) pp1))))
  344.       ((eq (var-l pp1) (var-l pp2))
  345.         (list (mkterm (var-l pp1) (max (exp-l pp1) (exp-l pp2))
  346.               (power-product-lcm (coef-l pp1) (coef-l pp2)))))   
  347.       ((greater-lead-p (vex-l pp1) (vex-l pp2))
  348.         (list (mkterm (var-l pp1) (exp-l pp1)
  349.               (power-product-lcm (coef-l pp1) pp2))))
  350.       ((greater-lead-p (vex-l pp2) (vex-l pp1))
  351.         (list (mkterm (var-l pp2) (exp-l pp2)
  352.               (power-product-lcm (coef-l pp2) pp1))))))
  353.  
  354.   (defun s-pol (p1 p2)
  355.     (let* ((lpp1 (leading-power-product p1))
  356.        (lpp2 (leading-power-product p2))
  357.        (lpc1 (leading-power-product-coef p1))
  358.        (lpc2 (leading-power-product-coef p2))
  359.        (lcpp (power-product-lcm lpp1 lpp2)))
  360.       (sub-pol
  361.         (multi-pol
  362.       lpc1
  363.       (multi-pol
  364.         (simple-power-product-quotient lcpp lpp1) p1))
  365.     (multi-pol
  366.       lpc2
  367.       (multi-pol
  368.         (simple-power-product-quotient lcpp lpp2) p2)))))
  369.         
  370.   (setq p1 (leading-power-product r1))
  371.   (setq p2 (leading-power-product r2))
  372.  
  373.   (export s-pol leading-power-product power-product-very-simply-divisible-p
  374.       power-product-lcm simple-power-product-quotient
  375.       leading-power-product-with-coef)
  376.       
  377. )
  378.  
  379.